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s Abstract 

9 In this paper we propose a procedure to evaluate Bayesian confidence 

10 intervals in counting experiments where both signal and background 

11 fluctuations are described by the Poisson statistics. The results obtained 

12 when the method is applied to the calculation of upper limits will also be 

13 illustrated. 

14 Keywords: Counting experiments, Confidence intervals, Upper limits 

15 1. Introduction 

16 The evaluation of confidence intervals and limits is a common task in 

17 particle physics [lj. Usually the goal of an experiment is that of determining 
ib a parameter 9 starting from a set of measurements of a random variable x (the 

19 outcome of the experiment) and assuming an hypothesis for the probability 

20 distribution function (p.d.f.) f(x\6). A confidence interval for 9 at the 

21 confidence level (or coverage probability) 1 — /3 includes the true value of 

22 9 with a probability 1 — (3. This means that if the experiment were repeated 

23 many times, the estimates of 9 will fall in the confidence interval in a fraction 
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24 1 — j3 of the experiments 0. 

25 Confidence intervals can be evaluated following either the frequentist 0, 

26 0, 0] or the Bayesian approaches (a detailed discussion about these methods 

27 can be found in ref . |l| ) . One of the main issues arising when applying the 

28 frequentist approach is that the confidence intervals may include unphysical 

29 regions for the parameter. These cases can be fixed by introducing some 

30 "ad-hoc" corrections in the mathematical procedures used to evaluate the 

31 intervals (see for instance ref 0,0]). Such issues can be easily avoided when 

32 following the Bayesian approach, by properly choosing the prior p. d.f. @ 

33 for the parameter. However, the fact that the Bayesian approach requires 

34 a choice of the prior p. d.f. for the parameter introduces some degree of 

35 arbitrariness in the evaluation of confidence intervals. 

36 In many cases the outcome of an experiment can be described in terms 

37 of a set parameters, not all being of any interest for the final result (nuisance 

38 parameters). In these cases the experimenter wishes to evaluate confidence 

39 regions on the parameters of interest, in a manner that is independent on 

40 the nuisance parameters. From a mathematical point of view, the Bayesian 

41 approach allows one to treat the nuisance parameters in a very simple and 

42 straightforward way. Indicating respectively with 9 and v the parameters of 

43 interest and the nuisance parameters, one has to write down their joint prior 

44 p. d.f. n(6,v)E and to evaluate from it the marginal prior p. d.f. for 9. In 

45 most cases 9 and v are independent random variables, and their joint p. d.f. 

46 can be factorized as tt(9, v) = 7i(9)tt(u), thus making the calculation easier. 

47 In the following sections we will illustrate an application of the Bayesian 

48 approach to the analysis of a counting experiment with background. We 

49 will assume that the outcome of the experiment can be modeled in terms of 

50 a parameter of interest (signal) and of a nuisance parameter (background), 

51 that will be supposed to be independent on each other. The posterior p. d.f. 



1 In the most general case 9 is a vector of parameters, i.e. 9 — (6i, 02, ■ • • , n ) and 
x represents a set of observables, i.e. x — (xi, x%, ■ • • , x m ). Consequently, the p. d.f. 
f(x\6) — f(xi, . . . , x m \9i, . . . , 9 n ) will be a function of m+n variables. Also, the confidence 
interval defined in the case of a single parameter becomes a n-dimensional confidence region 
for the vector of parameters 0. 

2 This can be done by assigning null probabilities to the values of 9 in the unphysical 
regions. 

3 Hereafter we will us the greek symbol 7r to indicate prior p.d.f.s and the latin symbol 
p to indicate posterior p.d.f.s 
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52 for the signal will be derived starting from a set of simple and quite general 

53 assumptions on the prior p.d.f.s for both the signal and the background. 

54 The formulas to evaluate upper limits on the signal will also be derived 

55 and the results will be discussed. Our results are a generalization of the 

56 ones illustrated in ref. [lj], where the same problem is discussed, but without 

57 taking background fluctuations into account. 

58 2. Formulation of the problem 

59 Let us consider a counting experiment in which one wants to measure the 
eo signal counts in presence of a background, that is measured in a subsidiary 

61 experiment. Such a situation can happen when measuring the activity of a 

62 weak radioactive source in presence of background with a counting device, 

63 like a Geiger-Muller detector. In this case two independent measurements 

64 are carried out, one for the signal and one for the background, and the results 

65 have to be properly combined. Another example is an experiment in which 
ee the signal and the background are evaluated looking at two different space 

67 regions. Also in this case, the signal is evaluated by combining the counts 

68 measured in the signal and background regions 0. 

69 In the following we shall denote with n and m the number of counts 

70 measured respectively in the signal and in the background regions. We will 

71 also indicate with s and b the true values of the signal and of the background 

72 counts respectively, and we will assume that the true value of the background 

73 counts in the signal region is given by cb, where c is a constant value that 

74 is assumed to be exactly known. Let us consider, as a first example, an 

75 experiment to measure the activity of a radioactive source. In this case 

76 n and m are respectively the counts recorded during a time interval T s in 

77 presence of the source, and during a time interval T& when the source has 

78 been removed. The value of c can be estimated as c = T s /Tb, and the source 

79 activity can be evaluated as s/T s once the value of s has been measured, 
so A second example, taken from the astrophysics, is the measurement of the 
si photon flux of a point source in the sky with a gamma-ray detector. In this 
82 case the signal region can be a cone centered on the source direction with a 



4 We prefer to use the word "combine" instead of "subtract" because in many cases 
the signal cannot be evaluated by simply subtracting the counts in the background region 
from the ones in the signal region. This may happen for instance when the counts in the 
background region are more than the ones in the signal region. 
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83 given angular aperture, while the background region can be an annulus far 

84 away from the source. In this case c can be evaluated as the ratio between 

85 the solid angles of the signal region and of the background region, eventually 

86 multiplied for the ratio between the live times of the two regions. 

87 Under the above assumptions, the probability of measuring m counts in 
as the background region will be a Poisson distribution with mean value b, i.e.: 

b m 

p(m\b) = e~ b — (1) 
mi 

89 while the probability of measuring n counts in the signal region will be a 

90 Poisson distribution with mean value s + cb, i.e.: 

p(n\ Sl b)=e-^ {S + ^ (2) 

91 Since the two measurements are independent, the joint p.d.f. for n and m 



92 



100 



102 
103 
104 
105 



will be given by: 



( s _j_ c b) n 

p(n } mis, b) = e~ {s+cb) - -— e~ 6 — (3) 

n\ ml 

93 Our problem is that of evaluating a Bayesian confidence interval (or 

94 an upper limit) for the parameter s, independently on b. An analogous 

95 problem is discussed in the textbook |6| and in ref . [1] , where the background 

96 value is assumed to be exactly known. Our discussion will be therefore 

97 a generalization of refs. [H, @]. A possible solution of the problem taking 

98 background fluctuations into account is given in ref. 0]. However, in ref. 0, 

99 a gaussian p.d.f. is used to model the background, the results being valid 
in the case of large counts. A similar problem is also illustrated in ref. jsj, 

ioi where Bayesian confidence interval are evaluated for a Poisson signal with 
known background and with fluctuations on the detection efficiency. 

Another possible solution to our problem is also given in refs. j3,[5j], where 
a frequentist approach is followed with the application of the profile likelihood 
method. However, while the procedure described in refs. 0, 11] requires some 
we adjustments to handle the cases when n < cm or when either n = or 
107 m = 0, in our method the treatement of these cases is straightforward and 
los does not require any adjustment. In fact, the formulas that will be derived 
109 in sections [3] and H] are valid for all the values of n and m. 
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no 3. The Bayesian approach 

in The implementation of the Bayesian approach requires the "probabilistic 

112 inversion" of eq. [31 i.e. the evaluation of the conditional p.d.f. p(s,b\n,m) 

in starting from p(n, m\s,b) and applying the Bayes'theorem. 

H4 3. 1 . Choice of the priors 

us The application of the Bayes'theorem starts from the assumption of a 

lie prior p.d.f. for both the random variables s and b. In the following we will 

n7 assume that s and b are independent. 

us For the true background value b we will assume a uniform prior: 

us where tt > is a constant. The parameter b is only constrained to be 

120 non-negative. 

121 On the other hand, for the signal true value s we will assume a prior p.d.f. 

122 given by: 

/ ks~ a s>0 
^ = \ s < (5) 

123 with k > 0. Also in this case, the only constraint on the parameter s is that 

124 it must be non-negative. 

125 It is worth to point out at this stage that both ir(b) and ir(s) defined in 

126 eqs. |H and are improper priors and then they can lead to posterior p.d.f.s 

127 that are not normalizable. In particular, as it will be shown in sec. I3.3[ this 

128 happens when setting a > 1. Our calculations will therefore be valid only 

129 for a < 1. 

130 The uniform prior for the signal is obtained by setting a = and 

131 represents the natural choice when the experimenter does not have any model 

132 for the signal. On the other hand, the choice of a power-law prior with a > 

133 will reflect the experimenter's belief that small signal values are more likely 

134 than larger ones. It is also possible to choose negative values of a: this choice, 

135 that is rather uncommon, would favour larger signal values with respect to 

136 smaller ones 0. 



5 It's easy to show that Tr(b)db — oo and 7r(s)ds = oo for any value of a. 

6 Let si < S2 be two possible signal values and let us consider the intervals [s±, si + As] 
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137 3.2. Evaluation of the background posterior p. d.f. 

138 Applying the Bayes'theorem and using for 7r(6) the expression in eq. H]it 

139 is possible to obtain the following equation: 

= KHjMjj = P(m\b) 
J p(m\b)ir(b)db J °° p(m\b)db 

uo Finally, replacing p(m\b) with its expression given in eq. HJ it is 

141 straightforward to obtain the final result: 

b m 

P(b\m) = e- b — (7) 
ml 

142 Note that even though the expression of p(b\m) in eq. [7] is the same as 

143 that of p(m\b) in eq. [TJ their meanings are completely different. In fact, the 

144 random variable in eq. [T]is m, and the formula tells that m follows a Poisson 

145 distribution; on the other hand, the random variable in eq. [7J is b, and the 
U6 formula tells that b follows a Gamma distribution. 

147 Using eq. [7J one can also easily evaluate the average value of b and its 
us standard deviation. It is easy to show that (b) — m + 1 and o~\> = y/m~+T. 

149 3.3. Evaluation of the signal posterior p. d.f. 

150 The Bayes'theorem can be applied to get the joint posterior p. d.f. for 

151 both s and b: 

p(n,m\s,b)n(s)n(b) 

p{s, b\n, m) = (8) 

J as J abp{n,m\s,b)7r{s)7r{b) 

152 Replacing p(n,m\s,b) with its expression given in eq. |3] and 7r(6) and 

153 7r(s) with their expressions given in eqs. [1] and [2] respectively, eq. |8] can be 

154 rewritten as follows: 

&K m) = / °° ds Jo 00 d6 e— (c+m™s-«(s + cbY (9) 



and [s2,S2 + As]. The ratio between the probabilities of finding s in the two intervals 
is given by R — P(si < s < si + As)/P(s 2 < s < s 2 + As) = (s 2 /si) Q . When a > 
(a < 0) then R > 1 (R < 1) and small (large) signal values are more likely than larger 
(smaller) values. For a discussion about the choice of priors in the Bayesian approach see 
for instance the textbook (a. 
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155 Indicating with N the denominator in the right-hand side of eq. [HI it can 

156 be rewritten as: 



/*oo roc /*oo 

N= / dbb m e^ c+1)b ds e- s s- a (s + cb) n = / db b m e~ {c+1)b f{b) (10) 
Jo Jo Jo 

157 where we have indicated with f(b) the result of the integral in ds, that can 

158 be seen as a function of the variable b. 

159 Applying the binomial theorem, the term (s + cb) n in the expression of 
wo f{b) can be expanded as follows: 



n | 



k=o K ' 



i6i Using this result, the expression of f(b) becomes: 



n (cb) n ~ k f°° 



k\(n-k)\ J 
k=o K ' Ji> 

162 and, taking into account the definition of the Gamma functional the previous 

163 equation can be written as: 

, F(k — a + 1) 

fib) = r n + 1) E cb mT7 W? ^TTT ( 13 ) 

^ r(fc + i)r(n - fc + 1) 

164 Introducing the expression of f(b) given by eq. [13] in the expression of iV 

165 given by eq. [TO] we get the following result: 



iv = r(n + i)^ 1 L— p- / d& 6 -+n- fee - (c+ i)6 14 
^ r(fc + i)r(n - fc + 1) Jo 



lee By making a proper change of variable, the integral in the right-hand side 

167 eq. can be expressed in terms of a Gamma function. Hence eq. O can be 

168 rewritten as follows: 



7 The Gamma function is defined as T(z) — J °° dte t t z 1 . It can be shown that 
r(z) = (z — 1)! when 2 is a positive integer. 
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r(n + l) ^ T(k-a + l)T(m + n-k + l) ( c Y~ 
( c + 1 ) m+1 t^ T(k + l)V(n-k + l) 

169 The joint posterior p.d.f. for s and b is then given by: 

p{s, b\n, m) = ^ e - s - {c+1)b b m S - a {s + cb) n (16) 

170 with the expression of N given in eq. [15] 

171 To evaluate the marginal p.d.f. for s we need to integrate the joint p.d.f. 

172 with respect to b: 



p{s\n, m) 



poo J pOO 

/ p(s,b\n,m)db= —e- s s~ a db e - (c+1)b b m (s + cb) n (17) 
Jo N Jo 



173 Indicating with g(s) the integral in the right-hand side of eq. [T71 it can 
be evaluated in a similar way to that used to calculate f(b). It is easy to 
show that: 



174 
175 



r(n + l) " T(m + n-k + l) ( c V " k 
9{S) ~ (c+l) m + 1 ^r(A; + l)r(r2-A; + l) \c+l) 1 j 

ne Introducing in eq. [T7] the expression of g(s) given by eq. [18] and the 

177 expression of iV given by eq. [151 h is possible to show that the posterior 

178 p.d.f. for s is given by: 

A r(m + v2-/c + l) / c \ sfc - Qe - s 

fc t^o T(A; + l)r(n - fc + l) 1 c+l) 

p{s\n,m) = (19) 

™ T(k - a + l)r(m + n - k + 1) / c 

fc ^ r(fc + l)r(n - k + 1) \cTl y 

179 It is worth to point out here that eq. [19] is valid only for a < 1. In fact, 
wo if a > 1, the argument k — a + 1 of the Gamma function in the sum in the 

181 right-hand side of the denominator will be null or negative, and consequently 

182 the Gamma function either will assume negative values or will diverge. As 

183 discussed in sec. 13.11 this behaviour is due to the fact that since the posterior 



8 



184 p.d.f. p(s\n,m) is obtained starting from improper priors, its normalization 

185 is not guaranteed. 

186 Fig. [I] shows the posterior p.d.f.s for the signal s evaluated for some 

187 different values of n and m in the case c = 1, i.e. when the background 
las region has the same size as the signal region. The calculations have been 
189 performed in the cases a = 0.5 (small signal expected), a = (uniform 
wo prior) and a = —0.5 (large signal expected). The differences between the 

191 three p.d.f.s become negligible when when n is larger than cm. 

192 From eq. [19] one can also easily calculate the moments of the p.d.f. 

193 p(s\n, m). It's easy to show that the hth moment of the p.d.f. is given 

194 by: 



202 
203 



, x n—k 

™ T(k - a + 1 + h)T(m + n - k + 1) 

ho r(fc + l)T(n - k + 1) \c + l, 

(s h \n, m) = V - ra 7 _ fc (20) 

n T(k - a + l)T(m +n-k+l) 

fc = r(Jfe + l)r(n - k + 1) V c + 1 

In particular, eq. [20] allows to calculate the expectation value of s, i.e. (s) 
and its variance, that can be evaluated as var(s) = (s 2 ) — (s) 2 . 

Finally, we can consider the limiting case with absence of background, 
that can be obtained by setting c = (background region with null size). In 
this limit, all the terms with [c/ (c + l)] n_fc in both the summations of eq. [T9l 
will vanish, but the ones with k — n, where \cj (c+ l)] n_fe = 1. Hence, in this 
case the posterior p.d.f. for the signal simplifies to: 

p(s\n,m) = — — (21) 

T(n — a + 1) 

and, as expected, does not depend on m. In particular, in the case a = 
the posterior p.d.f. for the signal becomes a Gamma distribution. 



204 4. Evaluation of the upper limits on the signal 

205 The posterior p.d.f. for the signal given by eq. [TH]can be used to evaluate 

206 Bayesian confidence intervals for s. In particular, in this section, we will 

207 apply the result of eq. [19] to the calculation of upper limits on s. 
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To evaluate the upper limit s u at the confidence level 1 — f3 we have to 

209 solve the integral equation: 

I p(s\n,m)ds = 1 -fi (22) 
Jo 

210 Taking advantage of the fact that 

s k ~ a e- s ds = 7(fc-a + l,s u ) (23) 



220 
221 
222 





211 where we have indicated with 7 the incomplete Gamma function 0, eq. 

212 can be rewritten as: 

/ \ n—k 

n j(k — a + 1, s u )T(m + n — k + 1) / 
ho T(k + l)T(n-k + l) I 7+ 1 

1_/3 = ; rr (2i) 

» T(k - a+l)T(m + n- k + 1) ( c 
ho T(k + l)T(n-k + l) \7+ 1 

213 Eq. EH can be solved numerically and allows to obtain the Bayesian upper 

214 limits on s at the confidence level 1 — (5 for any values of n and m, given 

215 the values of c and a. We have performed our calculations using the CERN 

216 ROOT package [9|. In particular, we have implemented a code that allows 

217 to evaluate the numerical solutions s u of eq. |24]for any given value of /3 with 

218 the Brent's method, using the ROOT built-in tools. 

219 In fig. [2] the upper limits on the signal at 90% confidence level are shown 
as a function of the observed counts in the signal (n) and background (m) 
regions in the case c = 1 for three different values of a. As expected, the 
choice of the power law index a in the signal prior p.d.f. will affect the result 

223 on the upper limits. In particular, the upper limits on the signal will increase 

224 with decreasing a. This result is in agreement with the fact that when a is 

225 positive and close to 1 small signal values are expected while, on the other 

226 hand, when a is negative, large signal values are expected. As pointed out 

227 in sec. (U this is a general issue of the Bayesian approach, the results being 

228 influenced by the initial belief of the experimenter. 



8 The incomplete Gamma function is denned as j(z,x) = f Q dte l t z - 1 . According to 
this definition 7(2, 00) = T(z). 
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229 Fig. E] shows the upper limits on the signal at 90% confidence level as a 

230 function of the background events m for several values of the signal events n. 

231 The calculation has been performed in the case c = 1 with the uniform prior 

232 (a = 0). For any given value of n, the upper limit on the signal decreases with 

233 increasing m, in agreement with the fact that (in the cases when n > cm) 

234 a rough estimate of s is given by n — cm and the upper limit is expected to 

235 be proportional to n — cm. It has also to be pointed out that in the case 

236 n = 0, i.e. when no events are observed in the signal region, the upper limit 

237 is always equal to 2.30, independently on m%. 

238 Finally, fig. 0] shows the upper limits on the signal at 90% confidence 

239 level as a function of the signal events m for several values of the background 

240 events m. Again, the calculation has been performed in the case c = 1 with 

241 the uniform prior (a = 0). The upper limits increase with increasing n, and 

242 the trend of the curves becomes almost linear when n > cm. 

243 Jf..l. Study of the frequentist coverage 

244 To study the frequentist coverage of the upper limits obtained with our 

245 procedure we implemented a dedicated simulation. For simplicity we studied 

246 only the case with c = 1, when the signal and the background regions have 

247 the same sizes. 

248 For any given pair of values of s and b, a sample of 10 5 events was 

249 simulated, each corresponding to the outcome of an experiment. Each 

250 event consists of a couple of random integer numbers (n, m), representing 

251 respectively the counts in the signal and in the background regions, that are 

252 generated according to the p.d.f. of eq. |3j For each couple of counts (n,m) 

253 the corresponding upper limits on s at 90% confidence level were evaluated 

254 by solving eq. [21] for different values of a. The coverage was then evaluated 

255 as the fraction of events with an associated upper limits less than the true 

256 value of the signal s. 

257 Fig. [5] shows, as an example, the results obtained with a simulated sample 

258 of events with s = 3.5 and 6 = 2. The distributions of the upper limits at 90% 

259 confidence level are shown for a = 0.5, a = and a = —0.5. As expected, 



9 It is evident that when no events (or a few events) are observed in the signal region 
and a large number of events are observed in the background region the prior p.d.f. given 
by eq. [5] is not adequate. In these cases a different prior p.d.f. should be used, that allows 
for negative values of s. 
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260 the coverage increases with decreasing a, since lower values of a correspond 

261 to more conservative upper limits. 

262 We have studied the frequentist coverage of our upper limits as a function 

263 of both s and 6. The results are summarized in Fig. [61 where the coverage is 

264 plotted as a function of s in the cases when 6 = 0,6 = 1,6 = 1.5,6 = 2,6 = 5 

265 and 6 = 10. As expected, the frequentist coverage decreases with increasing 

266 a. The choice of a uniform prior for the signal, i.e. a = guarantees a 

267 coverage that is larger than 90% for low signal values, and oscillates around 

268 90% with increasing s. On the other hand, the choice of a less conservative 

269 prior with a = 0.5 does not seems to affect significantly the coverage in case 

270 of an high background level. 



271 5. Conclusions 

272 We have developed a procedure that, following the Bayesian approach, 

273 allows to evaluate confidence intervals on the signal in experiments with 

274 background where both signal and background are modeled by the Poisson 

275 statistics. The implementation of the method is quite simple from the 

276 mathematical point of view, and does not require any adjustments to treat 

277 the cases when less events are observed in the signal region than those in 

278 the background region. The results obtained when our procedure is applied 

279 to the calculation of upper limits have been also discussed. The frequentist 

280 coverage of the upper limits evaluated with this procedure has been also 

281 studied. 
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Figure 1: The plots show the posterior p.d.f.s for the signal p(s\n,m) for different values 
of n,m and a in the case c = 1. Going from top to bottom, the plots correspond to 
n = 0,1,2,5; going from left to right the plots correspond to m — 0,1,2,5. The p.d.f.s 
have been evaluated assuming a = 0.5 (continuous lines), a = (dashed lines) and 
a = —0.5 (dotted lines). In each plot three vertical lines are also drawn with the same 
styles specified above, each indicating the average value of s from the corresponding p.d.f.s. 
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Figure 2: The plots show the values of the signal upper limits at 90% confidence level as a 
function of the counts observed in the signal and background regions, n and m, in the case 
c = 1. Going from top to bottom, the plots correspond to a = 0.5, a = and a = —0.5. 
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Figure 3: The plot shows the values of the signal upper limits at 90% confidence level 
as a function of the counts observed in the background region, to, in the case c = 1 and 
a = 0, for some different values of the counts in the signal region. The calculation has 
been performed for n = (•), 1 (■), 2 (a), 3 (▼), 5 (o), 10 (□) and 20 (A). 
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Figure 4: The plot shows the values of the signal upper limits at 90% confidence level as 
a function of the counts observed in the signal region, n, in the case c = 1 and a = 0, for 
some different values of the counts in the background region. The calculation has been 
performed for m = (•), 1 (■), 2 (a), 3 (▼), 5 (o), 10 (□) and 20 (A). 
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Figure 5: The plots show the distributions of the signal upper limits at 90% confidence 
level obtained from a sample of 10 5 simulated experiments with s = 3.5, 6 = 2 and c = 1. 
Going from top to bottom, the plots correspond to the upper limits evaluated by setting 
a = 0.5, a = and a — —0.5. The dashed lines indicate the true value of the signal. The 
coverage is graphically represented by the area of the histogram at the right of the dashed 
line normalized to the total area. 
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Figure 6: The plots show the frequentist coverage of the upper limits at 90% confidence 
level as a function of the true signal value s for different values of b with c = 1. From top 
left to bottom right the plots correspond to the cases b = 0, b = 1, b = 1.5, b = 2, 6 = 5 
and b = 10. The upper limits have been evaluated by setting a = 0.5 (continuous lines), 
a = (dashed lines) and a = —0.5 (dotted lines). The dash-dotted lines indicate the 90% 
coverage. 
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